#include <stdlib.h>
#include <cutil.h>
#include <stdio.h>
#include <cuda.h>
#include <shrUtils.h>

#include "util.h"
#include "matrixMul_kernel.cu"

#define BLOCK_SIZE 16

template <class T>
__host__ void initializeMatrix(T **, int, int);
template <class T>
__host__ void populateMatrix(T **, int, int);
template <class T>
__host__ void matrixMul(T**, T**, T**, int, int, int, const char*);
template <class T>
__host__ void deleteMatrix(T **, int);

int main(int argc, char** argv)
{
	//  Matrix dimensions
	int m = 0, n = 0, k = 0;

	m = atoi(argv[1]);
	n = atoi(argv[2]);
	k = atoi(argv[3]);
	
	
	int ** matrixA,   ** matrixB,   ** matrixC;
 
	matrixMul(matrixA, matrixB, matrixC, m,n,k, "int");
	
	//~ float ** matrixD,   ** matrixE,   ** matrixF;
	//~ 
	//~ matrixMul(matrixD, matrixE, matrixF, m,n,k, "float");
	//~ 
	//~ double ** matrixG,   ** matrixH,   ** matrixI;
	//~ 
	//~ matrixMul(matrixG, matrixH, matrixI, m,n,k, "double");	
	
	return 0;
	
}

//  Initializes matrix of size numRows x numCols
//  Values set to zero
template <class T>
__host__ void initializeMatrix(T **matrix, int numRows, int numCols) {
	for(int i = 0; i < numRows; i++)
		matrix[i] = new T[numCols];

	for(int i = 0; i < numRows; i++)
		for(int j = 0; j < numCols; j++)
			matrix[i][j] = (T) 0;
}

//  Populates matrix of size numRows x numCols
//  Values are floating point values between -10 and 10
template <class T>
__host__ void populateMatrix(T **matrix, int numRows, int numCols) {
	for(int i = 0; i < numRows; i++)
		for(int j = 0; j < numCols; j++)
			matrix[i][j] = (T) -10 + rand() * 20 / RAND_MAX;
}

//  Frees memory allocated to matrix with numRows rows
template <class T>
__host__ void deleteMatrix(T **matrix, int numRows)
{
    for(int i = 0; i < numRows; i++) {
            delete[] matrix[i];
            matrix[i] = NULL;
    }

    delete[] matrix;
    matrix = NULL;
}

template <class T>
__host__ void matrixMul(T** matrixA, T** matrixB, T** matrixC, int m, int n, int k, const char* datatype)
{
	//pointers to memory location on device
	T *  matrixA_p, *  matrixB_p, *  matrixC_p;
	int size = 0, device = 0;
	cudaDeviceProp deviceProp;
	
	unsigned int timer_execution = 0, timer_memcpy = 0;
	//  Allocate space for arrays
	matrixA = new T*[m];
	matrixB = new T*[n];
	matrixC = new T*[m];
	
	//printf("Initializing Matrices: A, B, C\n");
	initializeMatrix(matrixA, m,n);
	initializeMatrix(matrixB, n,k);
	initializeMatrix(matrixC, m,k);
	
	//printf("Populating Matrices: A, B, C\n");
	populateMatrix(matrixA, m,n);
	populateMatrix(matrixB, n,k);
	populateMatrix(matrixC, m,k);
	
	cudaSetDevice(1);
	
	cudaGetDevice(&device);         

	printf("device %d \n", device);
	
	cudaGetDeviceProperties(&deviceProp, device);
	printf("Using device %d: %s \n", device, deviceProp.name);

	
	//printf("allocating array on device\n");
	CUDA_SAFE_CALL( cudaMalloc((void**)&matrixA_p, sizeof(T)*m*n) );
	CUDA_SAFE_CALL( cudaMalloc((void**)&matrixB_p, sizeof(T)*n*k) );
	CUDA_SAFE_CALL( cudaMalloc((void**)&matrixC_p, sizeof(T)*m*k) );
	
	//printf("copying array to device\n");
	
	size = (sizeof(T)*m*n + sizeof(T)*n*k + sizeof(T)*m*k) / 1000000;
	
	printf("Moving %d MB: \n", size);
	
	CUT_SAFE_CALL( cutCreateTimer(&timer_memcpy) );
	CUT_SAFE_CALL( cutStartTimer(timer_memcpy) );
	

	CUDA_SAFE_CALL( cudaMemcpy(matrixA_p, matrixA, sizeof(T)*m*n, cudaMemcpyHostToDevice) );
	CUDA_SAFE_CALL( cudaMemcpy(matrixB_p, matrixB, sizeof(T)*n*k, cudaMemcpyHostToDevice) );
	CUDA_SAFE_CALL( cudaMemcpy(matrixC_p, matrixC, sizeof(T)*m*k, cudaMemcpyHostToDevice) );
	
	CUT_SAFE_CALL( cutStopTimer(timer_memcpy) );
	
	printf("time to transfer to device: %f seconds\n", cutGetTimerValue(timer_memcpy)/1000.0);
	printf("transfer bandwidth: %f MBps\n", (double)(size / (cutGetTimerValue(timer_memcpy)/1000.0)));
	
	//set the dimensions for the block and for the grid
	dim3 dimBlock( BLOCK_SIZE,     BLOCK_SIZE     );
	dim3 dimGrid(  k / dimBlock.x, m / dimBlock.y );
	
	printf("x: %d, y: %d\n", dimGrid.x, dimGrid.y);
	
	CUT_SAFE_CALL( cutCreateTimer(&timer_execution) );
	CUT_SAFE_CALL( cutStartTimer(timer_execution) );
	
	matrixMul_Kernel<<<dimGrid, dimBlock>>>(matrixA_p, matrixB_p, matrixC_p, m,n,k);
	cudaThreadSynchronize();
	
	CUT_SAFE_CALL( cutStopTimer(timer_execution) );
	printf("Execution time (%s): %f seconds\n", datatype, cutGetTimerValue(timer_execution)/1000.0);
	
	printf("freeing memory\n");
	
	deleteMatrix( matrixA, m );
	deleteMatrix( matrixB, n );
	deleteMatrix( matrixC, m );
	
	cudaFree( matrixA_p );
	cudaFree( matrixB_p );
	cudaFree( matrixC_p );
	
	CUT_SAFE_CALL( cutDeleteTimer(timer_execution) );
	CUT_SAFE_CALL( cutDeleteTimer(timer_memcpy) );
}
